library(Seurat)
library(tidyseurat)
library(harmony)
library(SingleR)
library(ggpubr)
library(tidyverse)
source('00_util_scripts/mod_seurat.R')

pathwang <- list.files('mission/HMCES_HIV/wang2020_GSE157829/', full.names = T) |>
  as_tibble() |>
  mutate(name = str_extract(value, 'C1|Q.'))

pathwang$name |> unique()

mtx3 <- c('matrix','barc','genes')

path3 <- mtx3 |>
  map(\(x)filter(pathwang, str_detect(value, x)) |>
        pull(value)) 

mtx_wang <- path3 |>
  pmap(\(x,y,z)ReadMtx(mtx = x,cells = y,features = z,strip.suffix = T),
       .progress = T)

mtx_wang |> glimpse()

# add idents to barcode
mtx_wang <- mtx_wang |>
  map2(unique(pathwang$name), add_name_field)

rownames(mtx_wang[[7]]) |> head()
rownames(mtx_wang[[7]]) <- rownames(mtx_wang[[7]]) |>
  str_remove('hg19_')

# unify gene name (rowname) before merging
united_genes <- mtx_wang |>
  map(rownames) |>
  reduce(intersect)

mtx_wang <- mtx_wang |>
  map(\(x)x[united_genes, ])

mtx_united <- mtx_wang |>
  reduce(RowMergeSparseMatrices)

sobj <- mtx_united |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj

sobj$mito.ratio <- sobj |>
  PercentageFeatureSet('^MT-')

sobj |> VlnPlot(c('mito.ratio','nFeature_RNA','nCount_RNA'),pt.size = 0)

sobj <- sobj |>
  filter(mito.ratio < 10)

sobj <- sobj |>
  mutate(batch = ifelse(orig.ident == 'Q7', 'hg19', 'hg38'),
               group = ifelse(orig.ident == 'C1', 'HC', 'HIV'))

# rds checkpoint ---------
sobj |> write_rds('mission/HMCES_HIV/wang2020hiv.rds')

sobj <- read_rds('mission/HMCES_HIV/wang2020hiv.rds')

sobj <- sobj |>
  quick_process_seurat()

hpca <- celldex::HumanPrimaryCellAtlasData()
monaco <- celldex::MonacoImmuneData()

sobj <- sobj |>
  mark_cell_type_singler(ref = hpca, new_label = 'hpca_main')

sobj |>
  DimPlot(cols = DiscretePalette(36), label = T, label.box = T)

sobj |>
  DimPlot(group.by = 'hpca_main')

sobj |>
  DotPlot(seurat_markers) +
  RotatedAxis()

sobj |>
  DotPlot('HMCES')

sobj |>
  DotPlot(cc.genes.updated.2019)
# 14 cluster may be cycling T cell?

sobj |>
  FindMarkers(ident.1 = 13, only.pos = T, logfc.threshold = 2) |>
  head()

GeneSymbolThesarus('IGJ')
# 13 cluster can be plasma cell / plasmablasts

sobj <- sobj |>
  mutate(manual_main = case_when(seurat_clusters == 13 ~ 'Plasma_cell',
                                 seurat_clusters %in% c(14,16) ~ 'Proliferating_T_cell',
                                 seurat_clusters %in% c(3,4,10) ~ 'CD8_T_cell',
                                 hpca_main == 'T_cells' ~ 'CD4_T_cell',
                                 .default = hpca_main))

sobj |> DimPlot(group.by = 'manual_main') +
  ggtitle('PBMC cell type')

sobj |> DotPlot('HMCES', group.by = 'manual_main')

# compare HIV vs HC in every leiden clusters
all_leiden <- sobj$seurat_clusters |>
  unique()

hmces_leiden <- all_leiden |>
  map(\(x)sobj |>
        FindMarkers(group.by = 'group', ident.1 = 'HIV',
                    features = 'HMCES',
                    subset.ident = x)) |>
  set_names(all_leiden) |>
  list_rbind(names_to = 'id')

hmces_leiden |>
  ggplot(aes(id, avg_log2FC, fill = p_val < .05)) +
  geom_col()

# compare HMCES in every cell type -------
all_type <- sobj$manual_main |>
  unique()

Idents(sobj) <- 'manual_main'

hmces_type <- all_type |>
  map(\(x)sobj |>
        FindMarkers(group.by = 'group', ident.1 = 'HIV',
                    features = 'HMCES',
                    subset.ident = x, logfc.threshold = 0), .progress = T) |>
  set_names(all_type) |>
  list_rbind(names_to = 'id')

hmces_logfc <- hmces_type |>
  mutate(manual_main = id, x = 1.5, y = 4,
         logfc = str_c('log2FC=',signif(avg_log2FC,2)))

man_meta <- sobj |>
  as_tibble() |>
  select(.cell, manual_main, group)

exp_hmces <- sobj |>
  get_abundance_sc_long(features = 'HMCES') |>
  left_join(man_meta)

exp_hmces |>
  ggplot(aes(group, .abundance_RNA, color = group)) +
  geom_violin() +
  geom_jitter(height = 0, width = .1) +
  geom_text(data = hmces_logfc, aes(x = x, y = 2.2, label = logfc), color = 'black') +
  scale_color_manual(values = c('blue','red')) +
  scale_y_continuous(expand = expansion(mult = c(0,0.2))) +
  theme_pubr() +
  facet_wrap(~manual_main, scales = 'free_y') +
  stat_compare_means(comparisons = list(c('HIV','HC')), label = 'p.signif', hide.ns = T) +
  labs(title = 'HMCES expression PBMC',
       subtitle = 'Wang 2020 GSE157829')


# compare HIV vs HC B cell
